Shell-Model Monte Carlo Simulations of Pairing in Few-Fermion Systems 
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We study a trapped system of fermions with an attractive zero-range two-body interaction using 
the Shell-Model Monte Carlo method. The method provides ab initio results in the low N limit 
where mean-field theory is not applicable. The energy and pairing properties are presented as 
functions of interaction strength, particle number, and temperature. In the interesting region where 
typical matrix elements of the two-body interaction are comparable to the level spacing of the trap 
we find large odd-even effects and signatures of shell structure. As a function of temperature, we 
observe the disappearance of these effects as in a phase transition. 
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The physics of ultracold atomic gases has been inten- 
sively pursued experimentally and theoretically in the 
last decade. Recently there has been great interest in 
strongly-interacting Fermi gases where Feshbach reso- 
nances allow the tuning of the two-body interaction. 
Studies of the transition from a dilute gas of fermionic 
atoms to a Bose condensate of molecules have therefore 
been possible in the laboratory Q, @, 0, 0, H, @]. While 
studies of degenerate Fermi gases have mostly dealt with 
large atom numbers and wide traps, efforts have begun 
to trap only few (1-100) atoms in tighter traps [7j. Also, 
with the implementation of three-dimensional optical lat- 
tices, a low-tunneling regime can be reached with es- 
sentially isolated harmonic oscillators containing only a 
few fermions at each site 0. This means that one can 
now explore few-body fermionic effects in trapped sys- 
tems with scattering lengths that are comparable to the 
inter-particle distance and the trap width. 

In this letter we report on a theoretical study of har- 
monically trapped fermions using the Shell-Model Monte 
Carlo (SMMC) approach. This method has been ex- 
tensively used in nuclear physics to determine nuclear 
properties at finite temperature in larger model spaces 
than can be handled in normal nuclear shell-model diag- 
onalization 0, [l(| ■ In the SMMC the many-body prob- 
lem is described by a canonical ensemble at temperature 
T = f3~ l and the Hubbard-Stratonovich transformation 
is used to linearize the imaginary-time many-body prop- 
agator e~@ H . Observables are then expressed as path 
integrals of one-body propagators in fluctuating auxil- 
iary fields. The method is in principle exact and subject 
only to statistical uncertainties. For equal mixtures of 
two hyperfine states at low density, the interaction can 
be modeled with an s-wave zero-range potential. It thus 
belongs to the class of two-body interactions which is free 



of sign problems in the SMMC [lOj . We present here the 
first application of this many-body method to ultracold 
gas physics. 

We study the usual zero-range pairing Hamiltonian 
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where we sum over all particles i and [ij] denotes a sum 
over fermion pairs with opposite internal states. The trap 
frequency is uj, and Vo denotes the interaction strength. 
The SMMC method was originally set up to handle nu- 
cleons, and pairing problems have been extensively con- 
sidered (ll| . The spin symmetric ultracold Fermi gas 
can now be easily mapped onto a single spin 1/2 nucleon 
species 12j, as the two-body interaction in the Hamilto- 



nian corresponds to standard nuclear pairing. 

Simple dimensional arguments reveal that the matrix 
elements scale as 1/fe 3 , where b — h/moj is the oscillator 
length. It is therefore natural to redefine the interac- 
tion strength in terms of Vo = —ghcob 3 , where g is a 
dimensionless strength measure. To relate Vo to the s- 
wave scattering length, one needs to regularize the zero- 
range potential and consider the finite model space cut- 
off, as pointed out in Here we will adopt the treat- 
ment in [l3 | and renormalize through continuum low- 
energy scattering. The energy cut-off used can be writ- 
ten E c = a 2 frui with a 2 — N max + 3/2. This defines a 
relation between Vo, a, and E c , which can be written as 
Ana/b = g/(ag/2ir 2 — 1). This gives the effective inter- 
action strength needed for a model space with N max + 1 
major shell. The relevant parameter regime is expected 
to be where the natural energy scale, given by the level 
spacing hu>, is comparable to typical two-body matrix 
elements between the trap states at T = (as a quanti- 
tative measure we use the (Is) 2 — (Is) 2 element). This 
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TABLE I: Energies (in units of TuuS) calculated with the 
SMMC method for a trapped fermion gas with interaction 
parameter g — 10 calculated at temperatures (in units of hu) 
T — 1/6 and T = 1/5 for different particle numbers N. The 
statistical uncertainty is given in parenthesis. HOSD denotes 
the non-interacting energies at T = 0. 
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HOSD T= 1/6 T= 1/5 
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2.60(1) 


2.62(1) 


12 


32 


29.8(2) 


29.9(1) 
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5.5 


5.12(1) 


5.14(1) 


13 


35.5 


32.8(4) 


33.3(2) 
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7.31(1) 


7.37(1) 


14 


39 


36.4(3) 


36.6(2) 
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10.5 


9.73(2) 


9.76(1) 
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42.5 


39.8(2) 


39.9(2) 
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13 


11.9(1) 


11.9(1) 


16 


46 


43.0(2) 


43.1(2) 
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15.5 


14.2(2) 


14.3(1) 
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49.5 


46.3(2) 


46.5(2) 
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18 


16.3(1) 


16.4(1) 
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53 


49.5(2) 


49.7(2) 
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21.5 


19.8(1) 


19.9(1) 
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56.5 
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53.1(2) 
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25 


23.1(2) 


23.2(1) 


20 


60 


56.1(2) 


56.4(1) 


11 


28.5 


26.5(2) 


26.6(1) 











turns out to be around g ~ 10 in our setup, which cor- 
responds to a = 10.666. The different regions of interest 
in terms of interaction strength and hu> in the context 
of atomic gas physics are discussed in 15, 1(| ■ We note, 
however, that our procedure to relate g and a is only 
approximate and modifications to the simple relation we 
use can occur, e.g. in tight traps [31 1. 

Previous works have considered few-fermion systems 
using advanced many-body methods. The Green's func- 
tion Monte Carlo methods were applied to homogeneous 
17| . as well as trapped systems 18j. No-core 19], and 
traditional shell- models [2fJ, using effective interactions 
have also recently been applied to these systems, particu- 
larly for very low particle numbers where exact results are 
available pll. [22I]. Finite- temperature, non-perturbative 
lattice methods have also been applied to homogeneous 
[23| and trapped fermions 24 1 . These works focus on the 
unitary \a\ — ► 00 limit and the crossover regime around 
it. The present SMMC approach is not variational, but 
an exact many-body method subject only to statistical 
uncertainties. Our choices of N max and g discussed above 
imply that the SMMC results given here are on the BEC 
side of the crossover regime. 

In Tab. U we present the SMMC energies for two dif- 
ferent temperatures small compared to hu>, along with 
the non-interacting energies. There is agreement between 
the two SMMC calculations within the uncertainties for 
almost all particle numbers. Notice that the absolute en- 
ergy uncertainty decreases with increasing T , as expected 
[lOj. The N = 2 problem can be solved exactly for all 
ajiil, and with a = 10.666 we find 2.4 for the T = 
ground-state. Our value of 2.6 at T = 1/&Tujj thus indi- 
cates that we are close to the T = limit. In general, 
we find that our energies are typically above those calcu- 
lated at unitarity in 3, 2l[ and in the crossover regime 
in 25j. This is expected due to finite temperatures and 




FIG. 1: (color online) The temperature dependence of the 
SMMC energies (both in units of hui) for selected even N 
systems and for strength parameters g — 10 (left) and 20 
(right). The zero-point energies are subtracted. Statistical 
uncertainties are too small to show. The dashed curves give 
the specific heat for N = 2, 8, and 10. Also shown are the 
exact N = 2 results (Busch et. al., [261 ]) with a = 10.666 
for the full space. The dash-dotted curve N = 2' gives the 
SMMC results for 5 major shells with g = 8.98. 



also because of Hilbert space differences (see discussion 
below) . 

Fig. Q] shows the SMMC energies as a function of tem- 
perature for selected even particle numbers and for the 
strength parameters g = 10 and 20. We can clearly see 
the convergence of our results at small T, approaching the 
T=0 ground state value. Quantitatively we find that the 
energies increase by about 1% (2%) going from T = Huj/6 
to T = hcu/5 (T = Tilj/A). At high T, the curves flat- 
ten out as they approach the equipartition limit of equal 
occupancies in the finite model space. 

We have used moderately small model spaces of 4 ma- 
jor shells for most calculations, which might be prob- 
lematic when using zero-range interactions. However, to 
check our convergence we have repeated the calculation 
for N=2 in an enlarged model space of 5 major shells with 
the interaction strength readjusted to the same scatter- 
ing length parameter. As can be seen in Fig. 1 the 
two calculations, performed for g = 10 for N max = 3 
and g = 8.98 for N max = 4, agree within the statisti- 
cal uncertainty up to a temperature T w 0.7. We also 
compare our SMMC results with the exact solution of 
the N=2 system as given in Ref. [25]. The comparison, 
although quite reasonable at low temperatures, can only 
be indicative as the pseudopotential used in [25] is simi- 
lar, but not identical to the present delta interaction and 
allows for a 2-body 1/r bound state which our calcula- 
tion does not include. At large T, the energy expectation 
value approaches E(T) — SNksT in the full model space. 
This behavior is, however, not obtained in our restricted 
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spaces where E(T) converges to the thermal equilibrium 
values which are 4.5 and 6 for N max = 3 and 4, respec- 
tively. 

The g — 10 results show a sharp increase in the energy 
around T ~ 0.3. Correspondingly the specific heat C — 
dE/dT exhibits a distinct peak around this temperature. 
Such a peak in the specific heat is usually associated with 
phase transitions in finite systems and indicates that the 
system changes its character. For g = 20 we see the 
peak at higher T as expected. As we will show below, it 
can be interpreted as the analogue of the superfluid-to- 
normal phase transition in finite nuclear systems 27, 28| 
and is associated with the breaking of pairs. Notice in 
particular that the peak height is larger for N = 8 than 
N = 10. This is a result of the shell closure at N = 
8 which causes the energy to be released in a smaller 
interval as T increases. The sizable width of the peak 
is due to the finite size of the system. The model space 
discussion above demonstrates that o ur p eaks are genuine 
and not a finite size Schottky effect [10( . 

To relate the peak in the specific heat to changes in the 
pairing correlations, we have calculated the expectation 
value of a number-conserving BCS-like pair matrix 



M, 



a, a' 



with the J = pair operator 



At = 



1 



x a. 



JM=Q0 
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where aj creates a particle in orbit j a (which is the com- 
bination of orbital and spin angular momentum of the 
fermions). This operator is thus a measure of the pairs 
with J — 0. An indication of the pairing correlations 
can be obtained from the sum over all matrix elements, 
called pairing strength in the following . Since we are 
only interested in genuine pair correlations we subtract 
the 'mean-field' values obtained at the same T but with 
.9 = 0. 

In Fig. [5] we show the pairing strength as a func- 
tion of particle number for different temperatures. The 
most striking feature is naturally the odd-even stagger- 
ing. The relative reduction of pairing strength for odd- 
particle numbers is related to the blocking of scattering 
of pairs into the orbital occupied by the unpaired parti- 
cle. The non-interacting systems have closed-shell con- 
figurations for N — 2,8, 20. With interaction switched 
on, these configurations manifest themselves by a relative 
reduction of the pairing strength (overlaid by a general 
increase due to a growing number of pairs) and a larger 
resistance against temperature increase. The strong dips 
observed for particle numbers N — 7, 9 and 19 are also re- 
lated to the shell closures. Relatedly the pairing strength 
is largest for mid-shell systems. From the inset we see 
that the staggering is larger for g — 20 and persists to 
larger temperatures as expected. 




FIG. 2: (color online) The pairing strength as a function of 
particle number, N, for various temperatures, T (in units of 
Tiui) for g = 10. The upper left inset shows the results for 
g = 20. The uncertainties are very small and not shown. 




FIG. 3: (color online) The pairing strength as a function of 
temperature, T (in units of hui), for various particle numbers, 
N and for g = 10 (left) and g — 20 (right). Uncertainties 
are small and not shown. Note the different scales in the two 
panels. 



As function of temperature, the pairing strength for 
g = 10 decreases and the odd-even effects vanish around 
T ~ (1/3 — X/2)Tkj except for the lowest N. This is con- 
sistent with the value at which the energies show sharp 
increases in Fig. [1] with resulting peaks in the specific 
heats. The same is seen for the g — 20 results. At higher 
T we see only a monotonous increase with N, indicat- 
ing equipartition in the model space with loss of all shell 
structure. 

To investigate further the transition between a paired 
state and a normal state, we show the pairing strength for 
g = 10 and g = 20 as a function of T for selected particle 
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numbers in Fig. [3] We note that the pairing correlations 
decrease rapidly for all particle numbers in the temper- 
ature regime which corresponds to the peak structure in 
the specific heat (with the strong decrease happening at 
larger T for g = 20 as expected) . We have calculated the 
pairing strength up to T = 4, at which point it has prac- 
tically vanished for all N. We believe this indicates that 
the transition temperature for our low particle numbers 
is some numerical factor less than unity (depending on g 
and N) times the trap spacing Tilo. Notice how the closed 
shell numbers N = 2 and 8 have different low-T behav- 
ior from the other systems (also the case for N = 20 
which is not shown) . Here the pairing remains relatively 
constant to larger temperatures due to the shell closure. 
Shown are also N = 7 and N = 9 which exhibit the large 
dips in the pairing strength in Fig. O They are gener- 
ally below the neighboring even-iV systems at low T, yet 
again confirming that the unpaired particle has a signifi- 
cant blocking effect on the pairing strength. Furthermore 
they have the same structure as the neighboring closed 
shell N = 8, but at lower magnitudes, ft is also interest- 
ing to observe that, for N = 7 and 9, the pairing strength 
is largest at finite T reflecting the competition between 
blocking by the unpaired particle and thermal excita- 
tions which moves the unpaired particle across the shell 
closure reducing the blocking effect. A similar effect has 
been found in SMMC studies for nuclei with odd nucleon 
numbers [3(| • At higher T these effects are smoothed out 
as pairing vanishes and the pairing strength orders with 
increasing number of particles. Comparing the g = 10 
and g = 20 results we see that the above effects are more 
pronounced for the stronger pairing strength and persists 
to higher temperature. This is consistent with the dis- 
cussions above. Similar evidence for a transition at finite 
T in the homogeneous s yste m in both energy and pair 



correlation was found in 24 1 



We have also considered the regime of g < 10 and 
found that at g = 1 (deep BCS regime) there are ba- 
sically none of the interesting effects left. We therefore 
predict that for small systems with simple two-body pair- 
ing, there is a regime at a > b where pairing features are 
pronounced. For g > 20 the calculation becomes unsta- 
ble, so we cannot access this deep BEC regime. 

In conclusion, the SMMC offers a good quantitative 
description of the behavior of equal mixtures of fermions 
in an interesting interaction regime. We studied the case 
of isotropic traps, and we extracted a number of par- 
ticularly relevant properties. The excellent convergence 
properties of the method holds promise for application 
to specific situations where, e.g., deformation properties 
and formation of higher angular momentum pairs may 
become relevant. Deformation was studied in the nu- 
clear case using the SMMC to predict shape transitions 
occurring as a function of temperature in the competition 
between pairing and quadrupole interactions |l0j . This 
is relevant also for ultracold gases with the recent realiza- 



tion of condensates with intrinsic long-range interactions 
between the atoms 32j. Recently the SMMC was also 
used to study the parity and spin properties of the den- 
sity of states in nuclear systems 3j| 34j. Transforming 
this to the atomic system could help us understand the 
low-energy excitation spectrum and the response of the 
gas to external perturbations. 

The authors wish to thank A.S. Jensen, M. Th0gersen, 
and A. Juodagalvis for many fruitful discussions. 
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